MINUIT_FIT

Overview

The MINUIT_FIT function performs nonlinear least-squares curve fitting using the MINUIT optimization algorithm, returning both fitted parameter values and their uncertainty estimates. This implementation leverages iminuit, a Python interface to the Minuit2 C++ library maintained by CERN’s ROOT team.

MINUIT (Function Minimization and Error Analysis) is a battle-hardened numerical optimization package originally developed at CERN for particle physics research. It was first described in the seminal 1975 paper by Fred James and Matts Roos and has been continuously refined over decades. The algorithm excels at minimizing statistical cost functions, particularly for maximum-likelihood and least-squares fits where reliable parameter uncertainty estimation is critical.

The function uses two core MINUIT algorithms:

MIGRAD: The primary minimizer, based on a variable metric method (a quasi-Newton approach using the Davidon-Fletcher-Powell formula). MIGRAD iteratively updates an approximation to the inverse Hessian matrix while searching for the minimum.
HESSE: Computes parameter uncertainties by calculating the full second-derivative (Hessian) matrix at the minimum and inverting it to obtain the covariance matrix.

For least-squares fitting, the function minimizes the chi-squared statistic:

\chi^2 = \sum_{i=1}^{n} \left( \frac{y_i - f(x_i; \boldsymbol{\theta})}{\sigma_i} \right)^2

where y_i are observed values, f(x_i; \boldsymbol{\theta}) is the model function with parameters \boldsymbol{\theta}, and \sigma_i are measurement uncertainties (optional). The reported standard errors correspond to the square roots of the diagonal elements of the covariance matrix at the minimum.

The function accepts arbitrary model expressions as strings, allowing flexible fitting to polynomial, exponential, trigonometric, and custom mathematical models. Parameter bounds can be specified to constrain the optimization search space. For more details, see the iminuit documentation and tutorials. The source code is available on GitHub.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=MINUIT_FIT(model_expr, xdata, ydata, param_names, param_guess, param_bounds, y_errors)

model_expr (str, required): The model expression as a string (e.g., “a * x[0] + b”). Use x[…] for predictors. Excel ‘^’ is converted to ’**’.
xdata (list[list], required): Predictor (independent) variables as a 2D list where each row is a sample (e.g., [[1.0],[2.0],…]).
ydata (list[list], required): Observed (dependent) values as a 2D list where each row is [[y_i]].
param_names (list[list], required): Parameter names used in the model expression as a 2D list (single row).
param_guess (list[list], required): Initial parameter guesses as a 2D list (single row).
param_bounds (list[list], optional, default: null): Optional parameter bounds. Either [[lower…],[upper…]] or per-parameter [[low,high],…]. Use None for unbounded.
y_errors (list[list], optional, default: null): Optional y uncertainties for weighted least squares as a 2D list [[err_1],[err_2],…]. If null, uses 1.0.

Returns (list[list] or str): 2D list [param_names, fitted_values, std_errors], or error message string.

Examples

**Example 1: Linear fit y = a*x + b**

Inputs:

model_expr	xdata	ydata	param_names		param_guess
a*x[0] + b	1	3	a	b	1	0
	2	5
	3	7
	4	9

Excel formula:

=MINUIT_FIT("a*x[0] + b", {1;2;3;4}, {3;5;7;9}, {"a","b"}, {1,0})

Expected output:

a	b
2	1
0.4472135955	1.2247448713

Example 2: Quadratic fit y = ax^2 + bx + c

Inputs:

model_expr	xdata	ydata	param_names			param_guess
ax[0]^2 + bx[0] + c	0	1	a	b	c	1	1	1
	1	2
	2	5
	3	10

Excel formula:

=MINUIT_FIT("a*x[0]^2 + b*x[0] + c", {0;1;2;3}, {1;2;5;10}, {"a","b","c"}, {1,1,1})

Expected output:

a	b	c
1	0	1
0.5	1.5652475841	0.9746794345

Example 3: Linear fit with parameter bounds

Inputs:

model_expr	xdata	ydata	param_names		param_guess		param_bounds
m*x[0] + c	0	0.5	m	c	1	0.5	0.5	0
	1	1.5					2	1
	2	2.5
	3	3.5

Excel formula:

=MINUIT_FIT("m*x[0] + c", {0;1;2;3}, {0.5;1.5;2.5;3.5}, {"m","c"}, {1,0.5}, {0.5,0;2,1})

Expected output:

"non-error"

Example 4: Weighted linear fit with y_errors

Inputs:

model_expr	xdata	ydata	param_names		param_guess		y_errors
a*x[0] + b	1	2	a	b	2	0	0.1
	2	4					0.1
	3	6					0.1
	4	8					0.1
	5	10					0.1

Excel formula:

=MINUIT_FIT("a*x[0] + b", {1;2;3;4;5}, {2;4;6;8;10}, {"a","b"}, {2,0}, {0.1;0.1;0.1;0.1;0.1})

Expected output:

a	b
2	0
0.03162277661	0.10488088484

Python Code

import math
import re
import numpy as np

try:
    from iminuit import Minuit
    from iminuit.cost import LeastSquares
except Exception:
    Minuit = None
    LeastSquares = None

def minuit_fit(model_expr, xdata, ydata, param_names, param_guess, param_bounds=None, y_errors=None):
    """
    Fit an arbitrary model expression to data using iminuit least-squares minimization with uncertainty estimates.

    See: https://iminuit.readthedocs.io/

    This example function is provided as-is without any representation of accuracy.

    Args:
        model_expr (str): The model expression as a string (e.g., "a * x[0] + b"). Use x[...] for predictors. Excel '^' is converted to '**'.
        xdata (list[list]): Predictor (independent) variables as a 2D list where each row is a sample (e.g., [[1.0],[2.0],...]).
        ydata (list[list]): Observed (dependent) values as a 2D list where each row is [[y_i]].
        param_names (list[list]): Parameter names used in the model expression as a 2D list (single row).
        param_guess (list[list]): Initial parameter guesses as a 2D list (single row).
        param_bounds (list[list], optional): Optional parameter bounds. Either [[lower...],[upper...]] or per-parameter [[low,high],...]. Use None for unbounded. Default is None.
        y_errors (list[list], optional): Optional y uncertainties for weighted least squares as a 2D list [[err_1],[err_2],...]. If null, uses 1.0. Default is None.

    Returns:
        list[list] or str: 2D list [param_names, fitted_values, std_errors], or error message string.
    """
    def _to_2d(x):
        return [[x]] if not isinstance(x, list) else x

    def _parse_col_vector(name, data, min_rows=2):
        data = _to_2d(data)
        if not isinstance(data, list) or len(data) < min_rows:
            raise ValueError(f"{name} must be a 2D list with at least {min_rows} rows")
        out = []
        for i, row in enumerate(data):
            if not isinstance(row, list) or len(row) == 0:
                raise ValueError(f"{name} row {i} must be a non-empty list")
            out.append(float(row[0]))
        return np.asarray(out, dtype=np.float64)

    def _parse_xdata(xdata):
        xdata = _to_2d(xdata)
        if not isinstance(xdata, list) or len(xdata) < 2:
            raise ValueError("xdata must be a 2D list with at least 2 rows")
        if not isinstance(xdata[0], list) or len(xdata[0]) == 0:
            raise ValueError("xdata rows must be non-empty lists")
        rows = []
        n_features = len(xdata[0])
        for i, row in enumerate(xdata):
            if not isinstance(row, list) or len(row) != n_features:
                raise ValueError("xdata must be a rectangular 2D list")
            rows.append([float(v) for v in row])
        return np.asarray(rows, dtype=np.float64)

    if Minuit is None or LeastSquares is None:
        return "Error: iminuit is not available in this environment."

    if not isinstance(model_expr, str) or model_expr.strip() == "":
        return "Invalid input: model_expr must be a non-empty string"

    # Convert Excel caret exponentiation to Python
    expr_str = re.sub(r"\^", "**", str(model_expr))

    try:
        x_arr = _parse_xdata(xdata)
        y_arr = _parse_col_vector("ydata", ydata)
    except Exception as exc:
        return f"Invalid input: {exc}"

    if x_arr.shape[0] != y_arr.shape[0]:
        return "Invalid input: xdata and ydata must have the same number of rows"

    # param_names
    param_names = _to_2d(param_names)
    if not isinstance(param_names, list) or len(param_names) == 0:
        return "Invalid input: param_names cannot be empty"
    pnames = param_names[0] if isinstance(param_names[0], list) else param_names
    if not pnames:
        return "Invalid input: param_names cannot be empty"
    pnames = [str(p) for p in pnames]
    n_params = len(pnames)

    # param_guess
    param_guess = _to_2d(param_guess)
    if not isinstance(param_guess, list) or len(param_guess) == 0:
        return "Invalid input: param_guess cannot be empty"
    p0_list = param_guess[0] if isinstance(param_guess[0], list) else param_guess
    if len(p0_list) != n_params:
        return f"Invalid input: param_guess must have {n_params} values"
    try:
        p0_vals = [float(v) for v in p0_list]
    except Exception as exc:
        return f"Invalid input: param_guess must be numeric ({exc})"

    # y_errors (optional)
    if y_errors is not None:
        try:
            yerr_arr = _parse_col_vector("y_errors", y_errors)
        except Exception as exc:
            return f"Invalid input: {exc}"
        if yerr_arr.shape[0] != y_arr.shape[0]:
            return "Invalid input: y_errors must match ydata length"
        if np.any(~np.isfinite(yerr_arr)) or np.any(yerr_arr <= 0):
            return "Invalid input: y_errors must be positive finite numbers"
    else:
        yerr_arr = np.ones_like(y_arr)

    # bounds (optional)
    lower = [-np.inf] * n_params
    upper = [np.inf] * n_params
    if param_bounds is not None:
        try:
            pb = _to_2d(param_bounds)
            # Format A: [[lower...],[upper...]]
            if isinstance(pb, list) and len(pb) == 2 and all(isinstance(r, list) for r in pb):
                if len(pb[0]) != n_params or len(pb[1]) != n_params:
                    return f"Invalid input: param_bounds must have {n_params} values for each bound"
                for i in range(n_params):
                    lb = pb[0][i]
                    ub = pb[1][i]
                    if lb is not None:
                        lower[i] = float(lb)
                    if ub is not None:
                        upper[i] = float(ub)
            # Format B: per-parameter [[low,high], ...]
            elif isinstance(pb, list) and len(pb) == n_params and all(isinstance(r, list) and len(r) == 2 for r in pb):
                for i in range(n_params):
                    lb, ub = pb[i]
                    if lb is not None:
                        lower[i] = float(lb)
                    if ub is not None:
                        upper[i] = float(ub)
            else:
                return "Invalid input: param_bounds must be [[lower...],[upper...]] or [[low,high],...]"

            for i in range(n_params):
                if not np.isfinite(lower[i]):
                    lower[i] = -np.inf
                if not np.isfinite(upper[i]):
                    upper[i] = np.inf
                if lower[i] > upper[i]:
                    return f"Invalid input: lower bound > upper bound for parameter {pnames[i]}"
        except Exception as exc:
            return f"Invalid input: param_bounds must contain numeric values ({exc})"

    # Safe eval environment (no builtins)
    safe_globals = {
        "np": np,
        "numpy": np,
        "math": math,
        "__builtins__": {},
    }
    safe_globals.update({
        name: getattr(math, name)
        for name in dir(math)
        if not name.startswith("_")
    })
    safe_globals.update({
        "sin": np.sin,
        "cos": np.cos,
        "tan": np.tan,
        "asin": np.arcsin,
        "arcsin": np.arcsin,
        "acos": np.arccos,
        "arccos": np.arccos,
        "atan": np.arctan,
        "arctan": np.arctan,
        "sinh": np.sinh,
        "cosh": np.cosh,
        "tanh": np.tanh,
        "exp": np.exp,
        "log": np.log,
        "ln": np.log,
        "log10": np.log10,
        "sqrt": np.sqrt,
        "abs": np.abs,
        "pow": np.power,
        "pi": math.pi,
        "e": math.e,
        "inf": math.inf,
        "nan": math.nan,
    })

    # Quick syntax check for expr
    try:
        compile(expr_str, "<model_expr>", "eval")
    except Exception as exc:
        return f"Invalid model expression: {exc}"

    def _model(x_matrix, *params):
        out = np.empty(x_matrix.shape[0], dtype=np.float64)
        for i in range(x_matrix.shape[0]):
            x_row = x_matrix[i]
            local_vars = {"x": x_row}
            for pname, pval in zip(pnames, params):
                local_vars[pname] = float(pval)
            y_pred = eval(expr_str, safe_globals, local_vars)
            y_val = float(y_pred)
            if not math.isfinite(y_val):
                raise ValueError("Model returned non-finite value")
            out[i] = y_val
        return out

    try:
        least_squares = LeastSquares(x_arr, y_arr, yerr_arr, _model)
        m = Minuit(least_squares, *p0_vals, name=pnames)
        m.errordef = 1.0

        # Apply limits
        for i, pname in enumerate(pnames):
            lb = None if lower[i] == -np.inf else float(lower[i])
            ub = None if upper[i] == np.inf else float(upper[i])
            if lb is not None or ub is not None:
                m.limits[pname] = (lb, ub)

        res = m.migrad()
        if not res.fmin.is_valid:
            return "Optimization failed: minimization did not converge"

        try:
            m.hesse()
        except Exception:
            pass

        values = [float(m.values[p]) for p in pnames]
        errors = []
        for p in pnames:
            try:
                err = m.errors[p]
                errors.append(float(err) if err is not None and math.isfinite(float(err)) else 0.0)
            except Exception:
                errors.append(0.0)

        # Normalize tiny values to 0.0 for stability
        values = [0.0 if abs(v) < 1e-12 else v for v in values]
        errors = [0.0 if abs(e) < 1e-12 else e for e in errors]

        return [pnames, values, errors]

    except Exception as exc:
        return f"Solver error: {exc}"

Overview

Excel Usage

Examples

Python Code

Online Calculator